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Abstract. We derive a variational formulation for thin film magnetization problems 
in type-II superconductors written in terms of two variables, the electric field and the 
magnetization function. A numerical method, based on this formulation, makes it 
possible to accurately compute all variables of interest, including the electric field, for 
any value of the power in the power law current-voltage relation characterizing the 
superconducting material. For high power values we obtain a good approximation to 
the critical state model solution. Numerical simulation results are presented for simply 
and multiply connected films, and also for an inhomogeneous film. 



Electric field formulation for thin film magnetization problems 



2 



1. Introduction 

A thin superconductor film in a perpendicular magnetic field is the configuration 
typical of experiments with superconducting materials, and is employed in various 
physical devices (SQUIDs, magnetic traps for cold atoms, etc.). Macroscopically, the 
magnetization of type-II superconductors is well described by eddy current models with 
critical state [H |2] or power law current-voltage relations. Solving these highly 
nonlinear eddy current problems helps to understand the peculiarities of magnetic flux 
penetration into thin films, and is necessary for the design of superconductor-based 
electronic devices. 

Analytically, the sheet current density is known for the Bean critical state model 
in both the thin disk [H [5] and strip [6] geometries. Numerical methods for modeling 
magnetization in flat films of arbitrary shapes were derived, for the power law model, 
by Brandt and co-workers in [3, E]; see also [9l [101 [H] and the references therein. For 
the critical state model a numerical scheme, based on a variational formulation of the 
thin film magnetization problems, has been proposed in [12]; see also [13]. Common 
to these numerical algorithms is the use of a scalar magnetization (stream) function as 
the main variable. The sheet current density in the film is obtained as the 2D curl of 
this function; the magnetic field can then be computed from the current density via the 
Biot-Savart law and compared to magneto-optical imaging results. 

The electric field in a superconductor is also of much interest: it is needed to find 
the distribution of the energy dissipation, which is often very nonuniform and can cause 
thermal instabilities. Computing the electric field e by means of existing numerical 
schemes can, however, be difficult for the power law model. 



where j is the sheet current density, j = Cq is a constant, jc is the critical sheet 
current density, and the power p is, typically, between 10 and 100. Indeed, even 
if the magnetization function is found with good accuracy, its numerical derivatives 
determining the sheet current density in the film are, inevitably, much less accurate. 
Computing the electric field via the constitutive relation increases the error further 
and makes it unacceptably large if the power p is high. 

As is well-known, the critical state model current-voltage relation can be regarded 
as the p oo limit of the power law (JT)); see [H] for a rigorous proof. The limit can be 
described as 



The electric field in this model can be nonzero only in a region where the current density 
is critical; there the field is parallel to current density and is determined by the eddy 
current problem with the constitutive relation ([2]). Note that even if the current density 
was computed, e.g., by means of the numerical scheme [12], the multi- valued relation 
(j2]) alone is not sufficient for the reconstruction of the electric field. Approximating the 




(1) 



IjI ^ jc] if \j\ < jc then e = 0; if e 7^ then e || j. 



(2) 
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electric field in a critical state model is relatively straightforward only for an infinite 
strip or a long superconducting cylinder in a perpendicular field [15] . 

For cylinders of an arbitrary cross-section in a parallel field, the magnetic field in 
the superconductor can be expressed via the distance to the boundary function (see 
[T7] ) or, in more complicated cases, found numerically. The current density is computed 
as the 2D curl of this field. Computing the electric field, however, remains non-trivial. 
A numerical algorithm for the electric field reconstruction, requiring integration along 
the paths of the magnetic flux penetration, has been proposed in [I6]. A dual/mixed 
variational formulation of magnetization problems served as a basis for the efficient 
computation of the electric field in [HI IT7] . 

Determining the electric field in thin film problems is more difficult. Under the 
simplifying assumption that the time derivative of the normal to the film magnetic 
field is, in the flux penetrated region, close to the ramping rate of the external field, 
approximate analytical expressions for the electric field were found for the Bean critical 
state model for the rectangular and related film shapes in [8]. 



Here we extend the approach [T8l [T7] and derive for thin film magnetization 
problems a convenient mixed variational formulation in terms of two variables: the 
electric field and a scalar auxiliary variable analogous to the magnetization function. 
We use Raviart-Thomas elements of the lowest order pTO] to approximate the (rotated) 
electric field and a continuous piecewise linear approximation for the auxiliary variable. 
Based on this approximation of the variational problem, our iterative numerical 
algorithm suffers no accuracy loss of the computed electric field even for very high 
values of the power p in ([1]). Hence, the algorithm can be used to find the electric and 
magnetic fields, and the current density for both the power and critical state model 
problems. 

In this work we focus on the derivation of the mixed variational formulation, 
describe the numerical algorithm, and present simulation results. Rigorous 
mathematical arguments, including the exact function space set up, and a proof of 
the algorithm convergence, etc., will be presented elsewhere |20] . 

2. Magnetization model: a mixed formulation 

Let Q G M."^ he a domain and, in the infinitely thin approximation, the superconducting 
film occupies the set {fi x 0} C M.^. By ei(xi,X2,t), where (xi,X2) G Q, we denote 
the tangential to the film (and continuous on it) component of the electric field 
e{xi,X2,Xs,t), and assume it is related to the film sheet current density j{xi,X2,t) 
by the power law ([1]). This law can be re- written as 



where r = 1/p and Ci = |ei|. The critical current density jc may depend only on 
{xi, X2) € Q (the Bean model for an inhomogeneous film) or also on the normal to the 
film component of the magnetic field (the Kim model). It is convenient to assume that 



j = 3c{e\/eoY ^Ci/eo 



(3) 
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Q is simply connected. If it contains holes, these can simply be filled in, with the sheet 
critical current density in the holes set to be zero or very small. 

In the outer space a; := \ x 0} we have the Faraday and Ampere laws, 

/^o + V X e = 0, V X h = 

with h\t=o — ho and h — )■ he(t) at infinity. Here /xq is the permeability of vacuum, e and 
h are the electric and magnetic fields, respectively, the given external magnetic field is 
uniform and normal to the film, he = (0, 0, hc{t)). We assume the initial magnetic field 
ho has zero divergence, V ■ /iq = 0, in u, its normal to the film component is continuous 
on {Q X 0}, and ho — he{0) = 0(|a;|"^) at infinity. 

We now relate the exterior space and film problems, then use the magnetic scalar 
potential and derive a 2D variational formulation, written for the electric field ei and 
the jump of magnetic potential on the film, which is convenient for the numerical 
approximation. 

The jump of the tangential component of the magnetic field across the cut {fl x 0} 
and the film current are related, 

j=n+x[h], (4) 

where n+ = (0,0,1). Here and below, [/] means the jump, — fla-, where 

= X {±0} are the two sides of {Q x 0}. Although the electric field e is not 
uniquely determined in oj by this model, its tangential component on the film, e-,-, is 
and has to be continuous: 

er|n+ = e^lf^- = ei. (5) 

Since in the outer space Vxh — 'Vx{h — he) — and fl is assumed to be simply 
connected, there exists a magnetic scalar potential w{x,t) such that h — he — — Vw. 
Furthermore, since V • — he) — V • h — 0, the scalar potential is a harmonic function 
in uj for any t, 

Aw = 0. (6) 

Integrating the Faraday law in time we obtain 

/xo (/i - /lo) + V X t/ = 0, (7) 

where U := f^edt' has the continuous tangential component Ur\n+ = Ur\n- = U\ : = 
ei dt' . Noting that the normal component of magnetic field is continuous on the film 
and also that n"*"- V x t/ = Curlt/i, where Curl / = dx^f2 — dx^fii obtain that 

lio n^-{h- ho) + Curl U; = 0. 

Substituting h = hc — Vw, we finalize our choice of the scalar potential w as the solution 
to the following exterior problem: 



Aw — in CO, 
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Si dt' 
w = 0{\x\-^) 



H 



on Q'^ and Q 



5 
(8) 



as X —J- cxD, 



with Ti = n+ ■ (he — Iiq). 

We set g = [w] and note that ii g = on the domain boundary dQ and is sufficiently 
regular (belongs to the space 5*0 = HqI'^{^1), see [HI I2U]), the unique solution to the 
following problem, 

Aw = in u, 
dw 
dn+ 

as Ixl —J- oo 



M = 9, 
w = 0{\x\''^) 



0, 



is the double layer potential ([2T], Ch. 3, §3.3 in the case of a closed surface Q, and [20] 
for the present choice of Q) 



w{x) 



1 

Air 



9{y) 



d 



dn: 



1 



F - y\ 



dy 



for x G cu, 



where d/driy = n+ ■ V^,. The normal derivative of this function, dw/dn~^, is continuous 
across the cut {Q x 0} and satisfies the variational equation 

dw 



5n+ 



ijj dx = —a{g, ip) 



(9) 



for any test function ip E Sq. Here the bilinear form 



1 

47r 
1 

47r 





n Jn 



Cur\g{x) ■ Curliply) 

\x - y\ 
Grad(y'(x) ■ Giadipiv) 
\x - y\ 



dx dy 



dx dy 



(10) 



is symmetric, and Curl0 = {dx^'t'i ~dxi(p) and Grad0 = {dxi(p,dx2(p) are 2D operators. 
We note that ^a{g, g) is the energy of the magnetic field induced by the ffim current 

J = ri+ X [h] = X [he — Vw] 
= Curl [w] = Curl (7. 

Substituting the normal derivative of w from ([8]) into the variational equation (Q we 
obtain 



'111 



for any ip G So] here {u,v)q, = J^uvdx is the inner product (or duality pairing) of two 
functions on Q. Differentiating with respect to time, we arrive at a more convenient 
form of this equation. 



a{dtg, ip)^ (Curl e;, ip)^ = - {dtK, ip)^ 



(12) 
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for any ip E Sq, with g\t=o = go determined by (ITT]) as 

aigo, i^) = -i (^c(O) - ^o) ■ n+, (13) 
for any ip E Sq. 

Finally, since j = Curl (7, we rewrite the current voltage relation ([3]) as 

Curl^ = jc(ei/eo)"~^ei/eo (14) 

and arrive at the mixed variational formulation ( fT2|) -([T4l) of the magnetization problem 
written for two variables, ei and g, defined on Q for any t > 0. 
It is convenient to use dimensionless variables, assuming 
_x' _ t' _ e[ 
L to Co 

. f , h' g' 

h = g 



jcO JcO jcoL 

where ' denotes dimensional physical quantities, 2L is the length of the projection of 
f2 onto the Xi-axis, the time scale to = fiojcoL/co, and jco is a characteristic value of 
the sheet critical current density. For homogeneous films with the field independent 
critical density jc we choose in our simulations jc/jco = 1- If this density depends on 
the normal to the film component of the magnetic field, jc = jc{hz) on {Vt x 0}, one can 
take ico = jc(0). The dimensionless form of the equations (|T2i) . is 

a{dtg, ip) + (Curl ei, ip)^ = - {dtK, (15) 

for any ip E Sq, and 

Curl^ = ^e[-^ei. (16) 

JcO 

Computing the normal to the film magnetic field component is needed in problems 
with field dependent critical sheet current densities, and also for the comparison of 
numerical simulation results to magneto-optical imaging. Noting that — = 
—dw/dri^ on {Vt x 0}, we can use ([9]) for determining the magnetic field component 
from the equation 

{hs- he,^)^ = a{g,ij) (17) 

for all ip E Sq. Alternatively, the explicit expression for dw/dri^ in (jS]) yields, in 
dimensionless variables, 

/i3 = /io3 - Curl f / Bi dt'\ . (18) 



Yet another possibility [I2] is to express the normal magnetic field component via the 
potential jump (magnetization function) g using the Biot-Savart law, 

h^{x,t) = K{t) + n+ ■ ^ / Vy(-^—A X j{y,t)dy 

47r Jo, \\x — y\j 

= K{t) - ^ I Grady (-^ — ^ — r ]■ Grady g{y,t) dy. (19) 
47r Jq \\x — y\/ 

These three approaches are further discussed in Sec. [31 
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3. Numerical scheme 

It is important to approximate the electric field ei in problem ( fT5l) - (|T6l) using curl 
conforming finite elements. In 2D problems, a simple change of variables leads to a 
formulation where curls are replaced by divergences; the divergence conforming Raviart- 
Thomas elements (see below) are an appropriate choice for such formulations. 
Let us substitute = Rv, where R is the rotation matrix 

1 
-1 

Taking into account that = |ej|, Curli? = — Div and i?^ Curl = Grad, we rewrite 
(II5])-([I6]) as 

aidtg, i^) - (Div V, = - {dth^, i^)^ , (20) 
Gradg = —\v\'-^v (21) 

jcO 

for any ip G 5*0. Here Divi; = dx^vi + 8x2^2 is the 2D divergence. Multiplying equation 
fl2T|) by a vector test function 77 and using Green's formula, we rewrite this equation as 

Jc"o'(jc|^^r~'«,r7)c + (^?,DivT7)f, = 0. (22) 
Equation (|T8l) should also be rewritten: 

/i3 = hoz + Div (^j^ V dt'^ . (23) 

We approximate f2 by a polygonal domain VL^. Let be a regular partitioning of 
QP' into triangles k and h = max^g-j-h diam(K) be their maximal size. Here vertices 
of lying on dQ^, the boundary of Q^, also lie on dfl. If Q contains subdomains 
with different critical current density values, the mesh is fitted in a similar way to the 
subdomain boundaries. By M'^ and S'^ we denote the sets of nodes and edges of this 
triangulation, respectively, with J\f^^ being the subset of the internal and J\f^ of the 
boundary nodes. Below, \X\ will denote the number of elements in the set X. 

Let 5^ be the space of continuous functions, linear on each triangle, and zero in the 
boundary nodes A/'^ • We define also the finite dimensional space of vectorial functions 
V^, linear on each triangle, v^\^ = + 6k(xi,X2), cl^ ^ I^^j G and such that the 
normal component of is continuous across any edge separating two adjacent triangles 
in T^. This is the space of divergence conforming Raviart-Thomas elements of the 
lowest order; see [19j for a detailed description of the edge related basis for V^. 

In addition, let = to < ^1 < • • • < ^at-i < ^at = T be a partitioning of [0, T] into 
possibly variable time steps t„ = t„ — t„_i, n = 1 ^ N. 
Our approximation of the problem (!20l) . ( l22l) is: 



Given G° e S^Joi n = l ^ N, find G and G such that 

a\G\^p'^)-r^ [Div V\^p'^)^,=a\G--\^P'^)-{h:-h:-\^p')^,, (24) 
JcO (Jc |V"r-i V^rj')' + (G",Divr7'^)^, = (25) 
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for all ip'^ G Sq and r]^ E . Here denotes he{tn), cl^{-,-) is defined by fllOj] with 
fl replaced by Q^, and {f,u)^ = YliK.eT'^^f averages the integrand / ■ u over each 
triangle k at its vertices: 

1 ^ 

m=l 

where {Pm}m=i the vertices of triangle n and \k\ its area. 

Furthermore, G Sq solves the corresponding approximation of (fT3!) . We note 
that it is not necessary to solve explicitly for as we can just replace the first term on 
the right-hand side of fl24l) for n = 1 by (/iqs — he{0),i>'^)nh. 

It is easy to show the existence and uniqueness of a solution to the nonlinear 
algebraic system fl21]) -f l2^ . see [SO]- To solve this system at each time level, we set 

= V"^^, denote \b\^ = \/\b\'^ + and approximate V" at the j^^ iteration 

by 

and find G"'^' G and V"-^' G such that 

= a'^(G"-\^^) - [K - K-^),t)^, (26) 

Jco' (Jc I V"'^ ry'^)' + (G"^ Divry'^)^, 

= Jc"o' (Jc ( I - ) , 77^^)" (27) 

for all ip^ G 5*9 and rj^ E . At each iteration, we need to solve the following linear 
system 

AG^ - TnBV^ = d, 
G^ + M^-^ = f-^ 

to determine G"'^ = Y.k=i^ij^h^GliJk and V"-'^ = Y.k=i^\£h\Vir}k; here {ipk} and {77^} 
are the standard bases for Sq and V^, respectively, and the time index n is omitted for 
simplicity. 

Here A is a symmetric positive definite full \J\fi^\ x \Afl^\ matrix with elements 
A.k,i = a^{'il)ki'4^i)'i M^~^ is a symmetric positive definite sparse \£^\ x \£^\ matrix with 
elements 

and 5 is a sparse \J^i^\ x \S^\ matrix with elements B^^i = {ipk,^^^ 'ni)n'^- 

We found that convergence of these iterations can be accelerated by supplementing 

them with an over-relaxation, i.e., by recalculating V^'^ as aV^'^ + (1 — a)V^'^~^ with 

a > 1 after each iteration. In all the examples below we chose e = 10^^ and a = 1.2. 
We also note that only the sparse matrix M^~^ must be recalculated at each 

iteration. The full matrix A is calculated only once for the chosen finite element mesh. 
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Since gradients of the basis functions ipk are constant on each triangle, to calculate these 
matrix elements one should find, see f lTOj) . the double surface integrals 

dx dy 




for every pair of triangles k;, G T^- We note that some of these integrals are singular. 
To accurately approximate this matrix, we followed the approach in the appendix of |22j : 
in particular, we used the exact analytical value [23] for the most singular cases / = m. 

To compare simulation results with the magneto-optical measurements of the 
normal to the film component of the magnetic field /is, an approximation H-^ to this 
component can be computed at the inner mesh nodes M^^ using a discretized form of 
the variational equation ( fTTll . 

IQ = K{tn)l + k'^ AG^, (29) 

where 1_ is the \N'^\ x 1 vector (1, 1, 1)'^ and A is the diagonal \Mi \ x \N'i \ matrix 
with Afc = ipk dx. Note that, due to the infinitely thin film approximation employed 
in our model, this magnetic field component becomes infinite on the domain boundary; 
see, e.g., the thin strip solution [6]. 

If the critical current density depends on the normal to the film magnetic field 
component, jc = jcihs), it is necessary to substitute jc by jc{H^'''~^) in fl28|) and to 
update the approximation after each iteration. The inner node values ( l29l) are 
not convenient for approximating in all triangles, as is needed in (128|1 . 

The piecewise constant approximation H^^ resulting from a discretization of ( 123|1 . 
can be written as 

M = H,, + C ijr-' + r„ ]/'■), (30) 

where i^oa is a piecewise constant approximation of /iqs, C is the sparse \T^\ x \£^\ 
matrix with elements Cm,i = Div(r7i)|K„, and [/"~^ denotes the coefficients of C/"~^ = 
Yl^Ji '^m V"^ for the standard basis of the Raviart-Thomas space V'^. We found, 
however, that such an approximation, based on the integrated over time electric field, 
can be inaccurate, especially in and around the film holes (in which the electric field 
remains undetermined in our model). 

A better piecewise constant approximation was obtained using a discretized form 
of f|T9|l : we set H_3\k '■= H^{pf^), where is the center of triangle k and 

Hiip.) = K - ^ E (GradG^) U, ■ / ds. (31) 

Here dn' is the boundary of k', i^^' is the unit outward normal to this boundary, and 
the integral over each side of triangle is computed numerically (as in we simply 
used Simpson's quadrature rule). 
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4. Simulation results 

The simulations have been performed in Matlab R2011a (64 bit) on a PC with Intel 
Core 15-2400 3.10Hz processor and 4Gb RAM. All film magnetization problems were 
solved for a growing external field he(t) = t and a zero field initial state. 

First, to test our method, we solved numerically the thin disk problem. Let Q be 
a circle of radius one. For the Bean critical state model the exact distribution of the 
sheet current density j = j{p,t) (p is known [H [5]. Here <p is the unit azimuthal vector 
in polar coordinates {p,(f))- In our dimensionless variables. 



where a{t) = 1/ cosh (2/;,e(t)) . The normal to the film component of the magnetic field 
can be found by means of a ID numerical integration using the equation 



where G{p,p') = K{k)/{p + p') - E{k)/{p - p'), k = 2^f^/{p + p') and K and E 

are complete elliptic integrals of the first and second kind. Furthermore, the electric 
field e = e{p,t)(f), where p^^dp{pe) = —dth^, e|p=o = 0. Approximating dths by 
{/t3(p, t„) — h2,{p,tn^i)}/Tn and integrating numerically, we calculate an approximation 
to the electric field distribution averaged over the time interval (tn-i,tn)- 

To compare with this semi-analytical solution of the Bean model, we set p = 1000 
in the power law model, and used our numerical algorithm flM|) -f l25|) with r = 1/p. Our 
numerical experiments confirmed that for a monotonically growing external field the 
current density and the magnetic field can be computed in one time step without any 
accuracy loss. The electric field e; is, however, determined by time derivatives of the 
external magnetic field and the magnetization function; in the discretized formulation 
f l2^ -f l25|) the field E"' = R V"' can be considered as an approximation to the field e; at 
some time moment in the interval or to the average field in this time interval. 

Hence, our numerical strategy was to make a large time step ti followed by a much 
smaller step T2 to obtain accurate approximations to all variables, including the electric 
field, at time t = ti + T2- The simulation results for t = 0.5 (Fig. [1]) were obtained 
with Ti = 0.45, T2 = 0.05 and two finite element meshes, with 4200 elements {h = 0.05) 
and 12000 elements {h = 0.03). Solution for these two meshes took, respectively, 5 
and 57 minutes of CPU time, not including the time for computing the full matrix A. 
Comparing to the solution of the Bean model described above, we found that the critical 
current zone was 0.65 < p < 1, where a(0.5) = 1/ coth(l) ^ 0.65, and the relative errors 
for the current density, the electric field, and the normal magnetic field component were, 
respectively, 1%, 4.9%, and 2.5% for the crude mesh and 0.6%, 3.1%, 1.4% for the fine 
mesh. Here the electric field in the semi-analytical solution for the Bean model was 
calculated using the hs distributions at the same two time moments, ti = 0.45, t2 = 0.5. 




} 



a{t) < P < 1, 
< p < a(t). 
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At each time level, the approximate current density was computed as CurlG", constant 
in each triangle. However, for the comparison we used the node values calculated, at 
each node, as the weighted-by-areas mean of the values in triangles to which the node 
belongs; such averaging increased the accuracy. We note also that the magnetic field 
/13 was determined using equation (l29l) . Hence, the field was found and compared to 
the exact solution at the internal nodes only (at the boundary nodes the exact field is 
infinite). 

In the next two examples we also assumed p = 1000, so the numerical solutions 
obtained should be close to solutions to the Bean model; the magnetic field was 
computed using the discretized Biot-Savart formula (1311) . 

The electric field is known to be strong near the film boundary indentations and, 
especially, in the vicinity of concave film corners (see Fig. |2]). Although similar problems 
have been solved by other authors before, this was done for p = 9 and p = 19 in |8] 
and [Si [10], respectively (and also for the critical state models in [12], but there without 
computing the electric field). Here we took time steps 0.3 + 0.08 + 0.02 = 0.4. 

In inhomogeneous films the electric field near the boundaries between regions of 
different critical current densities can be orders of magnitude higher than in other parts 
of the film (see Fig. [3]). Here jc/jco = 0.5 inside the rectangle, and jc/jco = 1 outside. 
We took time steps 0.2 + 0.2 + 0.1 = 0.5. The magnetic fiux penetrates deeper into the 
lower critical current density area. 

To solve a problem with a multiply connected film (Fig. ^ we filled the holes in 
and set jc/jco = 0.002 there, while keeping jc/jco = 1 in the film itself; we recall the 
electric field in the holes is not determined. This example was solved for p = 100 with 
time steps 0.47 + 0.03 = 0.5. A strong electric field is generated along the paths of fiux 
penetration into the holes. For p = 19 such problems were solved by a different method 
in [10]. 

5. Conclusion 

Existing numerical methods for thin film magnetization problems in type-II 
superconductivity are based on formulations written for one main variable: the 
magnetization function. The sheet current density, determined numerically as the curl 
of this function, is prone to numerical inaccuracy. The inaccuracy, usually tolerable in 
the current density itself, inhibits evaluation of the electric field by substituting this 
density into a power current-voltage relation if the power is high. For critical state 
models such an approach for computing the electric field is not applicable. 

The new variational formulation of thin film magnetization problems proposed in 
this work is written for two variables, the electric field and the magnetization function. 
The formulation serves as a basis for the approximation and computation of all variables 
of interest: the sheet current density and both the electric and magnetic fields. Our 
numerical algorithm remains accurate for any value of the power in the power law 
current- voltage relation. For high powers we obtain a good approximation to the solution 
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Figure 1. Thin disk in the perpendicular field, hc{t) = t. The Bean model solution 
(H [S] (black line) and the numerical solution (red dots) obtained with p ~ 1000, 
h = 0.03. Shown for t = 0.5: top - the modulus of the current density j; middle - the 
modulus of the electric field ei; bottom - the normal component of the magnetic field. 



of the Bean model. Evaluation of the local heat dissipation distribution in a film for 
both the power law and critical state models becomes straightforward. 

In this paper, we presented numerical simulation results for isotropic models 
with field independent critical sheet current density. However, our approach can be 
generalized to thin film problems with field-dependent [2] and anisotropic [23] sheet 
critical current densities. 
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Figure 3. Inhomogeneous film of elliptic shape in a growing external field; p = 1000. 
Sheet critical current density jc/ico = 0.5 in the rectangle and jc/jco = 1 outside of it. 
The finite element mesh contained 10,600 triangles and was refined near the boundary 
between the two regions (the blue line). Shown for t = 0.5: left - the modulus of the 
electric field ei; right - current streamlines (top) and levels of the normal to the film 
magnetic field component /13 (bottom). 

superconducting tapes Supercond. Sci. Technol. 24 075012 
[16] Badia-Majos A and Lopez C 2004 Electric field in hard superconductors with arbitrary cross 

section and general critical current law J. Appl. Phys. 95 8035-8040 
[17] Barrett J W and Prigozhin L 2010 A quasi-variational inequality problem in superconductivity 

Mathematical Models and Methods in Applied Sciences 20 679-706 
[18] Barrett J W and Prigozhin L 2006 Dual formulations in critical-state problems Interfaces and Free 

Boundaries 8 347-368 

[19] Bahriawati C and Carstensen C 2005 Three Matlab implementations of the lowest-order Raviart- 

Thomas MFEM with a posteriori error control Comput. Methods Appl. Math. 5 333-361 
[20] Barrett J W and Prigozhin L Existence and approximation of a mixed formulation for thin film 

magnetization problems (in preparation) 
[21] Nedelec J-C 2000 Acoustic and Electromagnetic Equations: Integral representations for Harmonic 

Problems ( Springer- Verlag: NY) 
[22] Sokolovsky V, Prigozhin L and Dikovsky V 2010 Meissner transport current in flat films of arbitrary 

shape and a magnetc trap for cold atoms Supercond. Sci. Technol. 23 065003 
[23] Arcioni P, Bressan M and Perregrini L 1997 On the evaluation of the double surface integrals 



15 




-1 1 



Figure 4. Circular film with three holes in a growing external field; p — 100. The 
finite element mesh contained 10,400 triangles and was refined near the domain and 
hole boundaries (blue lines). Shown for t — 0.5: left - the modulus of the electric field 
Ci; right - current streamlines (top) and levels of the normal to the film magnetic field 
component (bottom). 
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